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Abstract 

It  is  well  known  that  a  common  notion  of  time  in  distributed  systems  can  be  used  to  ensure 
additional  properties  such  as  real-time  behavior  or  the  identification  of  the  order  of  events.  As 
large-scale  hardware  testbeds  for  such  systems  are  neither  efficient  nor  easy  to  manage,  discrete 
event  simulations  (DES)  can  be  used  to  model  such  networks.  However,  to  ensure  an  exact 
behavior  of  such  simulations,  high  precise  models  of  the  local  clocks  are  also  needed:  the 
driving  oscillators  have  to  be  modeled  in  a  way  that  a  DES  simulation  of  a  free-running  node 
clock  shows  the  same  Allan  deviation  as  the  simulated  counterpart.  This  paper  shows  an 
approach  to  find  a  corresponding  model  for  a  simulator  using  white  noise  and  a  filter  with  the 
same  power  density  spectrum  as  real-world  oscillators. 


INTRODUCTION 

The  advantage  of  synchronized  clocks  for  network  nodes  is  very  well  known.  One  of  the  most  popular 
applications  is  time  division  multiple  access  (TDMA)  on  shared  media  and  the  direct  use  for  real-time 
tasks.  Besides  that,  the  simple  need  for  the  capability  of  sequential  ordering  of  events  occurring  in  large 
distributed  systems  [1]  depends  on  synchronized  clocks.  Standards  like  the  upcoming  IEEE  1588 
protocol  allow  the  synchronization  of  clocks  over  packet-oriented  networks. 

In  such  commercial  or  industrial  environments,  the  number  of  nodes  often  reaches  a  number  of  hundred 
or  more.  Consequently,  for  design  and  development  of  control  loops,  devices  or  even  enhancements  of 
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protocols,  the  implementation  of  prototype  testbeds  is  neither  cost  efficient  nor  trivial  to  manage.  Thus, 
there  is  a  need  for  simulation  of  the  operation  and  performance  of  such  large  systems. 

Besides  the  cost-efficient  implementation  of  a  simulation  for  large-scale  networks,  it  is  possible 
to  run  simulations  faster  or  slower  than  in  real  time1,  which  opens  new  possibilities  to  verify 
effects  that  are  too  difficult  to  see  in  hardware -based  testbeds.  Finally,  simulation  of  timely 
behavior  of  geographically  distributed  elements  has  the  advantage  of  having  a  common  notion  of 
real  time  available  at  each  and  every  node.  Such  fully  synchronized  clocks  at  every  node  can,  of 
course,  not  replace  the  actual  notion  of  time  at  a  node  (which  is,  of  course,  influenced  by  the 
local  clock  and,  therefore,  in  general  different  at  each  and  every  node). 

The  remainder  of  this  paper  is  structured  as  follows:  after  an  introduction  into  the  problems  of 
discrete  event  simulation,  the  theory  of  power  density  spectra  of  oscillators  will  be  given.  This 
will  be  followed  by  an  example  of  parameter  extraction  of  a  commercial  off-the-shelf  oscillator 
and  a  conclusion,  as  well  as  an  outlook  for  further  work. 

Discrete  Event  Simulators 

In  general,  two  possibilities  for  the  simulation  of  time-driven  systems  exist: 

•  Continuous  time  simulators,  where  a  timescale  is  cut  into  arbitrary  small  parts  and  the  status  of 
the  system  is  fully  determined  at  each  and  every  microtick  of  the  simulator.  However,  for  the  case 
of  high-precision  clock  synchronization,  this  turns  out  to  be  neither  efficient  nor  manageable  with 
today’s  reasonably  available  computing  power. 

•  Discrete  event  Simulation  (DES).  This  class  of  simulators  is  built  on  events,  which  can  be  raised 
at  a  certain,  so-called  simulation  time.  This  simulation  time  is  a  global  equally  spaced  timescale, 
such  as  TAI.  The  main  advantage  comes  from  the  fact  that  the  simulator  keeps  track  on  all 
upcoming  events  and  their  occurrence,  which  allows  the  system  to  “jump”  from  one  event  to  the 
next.  The  time  in  between  is  not  simulated  and,  therefore,  skipped  out,  which  increases  the 
performance  dramatically. 

When  using  DES  for  simulating  networks  performing  clock  synchronization,  each  simulated  node  has  to 
use  its  own  notion  of  time  and  not  the  globally  available  (and,  therefore,  evidentally  synchronized) 
simulation  time.  The  task  of  translating  this  simulation-wide  identical  time  notion  into  a  node-specific 
one  is  done  by  oscillator  models,  which  are  the  actual  topic  of  this  investigation. 

The  obviously  most  simplistic  approach  is  to  model  an  oscillator  with  a  drifting  frequency  offset. 
Although  this  model  has  shown  to  be  very  useful,  such  a  jitter-less  approach  turns  out  not  to  be  enough 
precise  for  simulation  scenarios  where  clock  synchronization  in  the  sub-nanosecond  range  is  needed.  Due 
to  the  fact  that  external  events  might  also  trigger  a  simulated  node  to  wake  up  and  take  action,  it  cannot  be 
assumed  that  there  is  no  earlier  run  (externally  inserted  new  event)  than  the  locally  calculated  next  event 
for  a  node. 

Problem  Statement 

This  non-continuous  progress  in  time  requires  special  handling  in  the  implementation  of  an  oscillator 


1  In  this  paper,  the  notation  real-time  shall  address  the  behavior  of  a  system  to  react  within  guaranteed 
(timely)  deadlines,  whereas  real  time  will  be  used  to  indicate  a  physical  timescale  such  as  TAI. 
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model.  This  means  that  any  oscillator  has  to  store  the  upcoming  future  of  his  scheduled  tick  events,  in 
order  to  produce  timely  consistent  output  in  case  of  an  event  insertion  (in  the  order  before  already 
scheduled  events). 

The  problem  with  the  simple  drifting  frequency  model  is,  besides  the  poor  mapping  of  physical  effects 
like  aging,  that  noise  and  statistical  effects  are  not  covered.  Consequently,  an  oscillator  model,  which 
gives  a  more  realistic  view  of  the  reality,  is  needed. 

While  the  determination  of  the  next  event  is  easy  when  using  a  non-statistical  behavior,  an  appropriate 
noise  generator  is  required  as  base  for  calculating  statistical  properties.  Common  simulators  provide 
random  number  generators  for  white  noise.  The  goal  for  the  design  of  the  presented  model  is  to  develop  a 
fdter  structure,  which  fdters  white  noise  in  a  way  that  the  Allan  deviation  of  real  and  modeled  oscillators 
match.  In  order  to  obtain  the  necessary  parameters  for  the  construction  of  a  matching  model,  the  power 
density  spectrum  and  the  setup  for  measuring  a  quartz  crystal  oscillator  is  used.  After  extracting  the 
required  parameters  (bandwidth,  border  frequencies,  and  the  comer  frequencies)  from  the  measurement, 
the  mathematical  construction  of  the  filter  completes  the  development  process. 


OSCILLATOR  MODEL 


The  Allan  variance  [2]  of  of  a  power  density  spectrum  Sy(j)  limited  to  \fic„fuc\  is  defined  by: 


°;X'pf,cJuc)-=  2J 


AsinJ^) 

A  {n fr)2  y 


(1) 


It  is  well  known  [3]  that  a  model  of  the  power  density  spectrum  of  an  oscillator  is  given  by: 


SVJ(f) 


_  k/7;c£/sA, 

(0  otherwise 


(2) 


with  Sy(f)=T*Syj(f)  with  -2  <  i  <  2.  Using  this  equation,  the  contributing  parts  of  the  Allan  variance  crv/ 
can  be  defined  as: 


^;ArJ,c'Juc)-=  2J 


sin4(^fr) 


S„i(f)df 


yd ' 


with  the  trivial  relationship: 

(r;  fic  =  X  ay,i  (T>  fic  ’fuc)- 

i=-2 


(3) 


(4) 


A  filter  with  a  pulse  response  that  satisfies  this  equation  will  respond  to  white  noise  with  an  output  having 
the  same  Allan  deviation  as  the  measured  system.  However,  the  evaluation  of  the  above  equations  for  fc 
— >  0  and  fuc  —>  °°  deliver  meaningful  results  only  for  values  of  i  =  {-2, -1,0},  whereas  the  parts  for  1  and  2 
reveal  infinite  results  for  the  integral. 

This  fact  is  well  known,  as  Jacques  Rutman  wrote  in  Characterization  of  Frequency  Stability  in  Precision 
Frequency  Sources  [4]:  “ It  is  necessary  to  specify  the  value  and  shape  of  the  low-pass  fdter  for  noise 
types  with  i  >  0.  Most  common  shapes  are  the  infinitely  sharp  low-pass  filter  ...” 
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Nevertheless,  for  values  of  fuc  <  oo,  a  solution  of  the  integrals  can  be  given  with  the  help  of  the  definition 
of  the  sine  integral: 

S(g)-.=  \l^u  (5) 


and  cosine  integral: 


C(f)  :=  y  +  tots-)  +  =  -f^du 


(6) 


[51,  as  well  as  the  Euler-Mascheroni  Constant  y  «  0.577  [61,  one  can  find  new  solutions  of  the 
contributory  parts  of  the  Allan  deviation  after  a  short  calculation. 


(J2y_2(r,0,fc)  =  7rr 


r8S(4iO-4S(2iO 


t 

sin2(ug)[4ug  +2ug  sin(2«g) + (8wg  -l)cos(2«g)  +  l]l 


3  u: 
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(7e) 
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&v,2  (r;0,fc)  = 


1  12w  +  sin(4w  )  -  8sin(2w  ) 


(kt) 


16 


(7b) 

(7c) 

(7d) 

(7e) 


with  the  substitution  «c(x)  =  nfc  r .  Finally,  the  actual  values  can  be  calculated  with 


crlM-fcfuc)  =  OyArAfJ)  -  o’LC^O./fc)  • 


(8) 


With  these  equations,  it  can  be  seen  that  the  parts  for  i  =  -2,-1,  and  0  are  independent  of/),  in  the  case  of 
ug  »  1.  Moreover,  for  i=l  and  2,  the  respective  parts  of  the  Allan  deviation  approach  infinity  with 
growing  numbers  of fc. 

Parameter  Extraction  from  Measured  Allan  Deviations 

The  parameter  extraction  now  can  be  done  very  easily  with  the  following  approach: 

1. )  The  values  of  a2y,i  and  a2y>2  are  limited  by  the  test  setup.  With  that  limitation,  the  lower  and  the 

upper  cut-off  frequency,  which  are  defined  by  the  measurement  equipment,  can  be  defined. 
Moreover,  the  relationship  2  nfuc  x  »  1  should  be  satisfied  all  relevant  values  of  x. 

2. )  The  graphical  curve  of  a2y  is  measured  empirically. 

3. )  The  theoretical  values  of  the  variances  (in  multiples  of  h,)  have  to  be  determined  via  equation  (7). 
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4.)  The  values  of  /?,  get  determined  numerically  with  the  method  of  the  infinity  norm  ||  .  ||^. 

With  that  parameters  the  power  density  spectrum  of  Sy(f)  is  sufficiently  determined,  and  can  be  used  for  a 
filter  design. 


OSCILLATOR  CHARACTERIZATION 

As  the  study  of  the  long-term  stability  of  an  oscillator  is  rather  time  consuming,  it  is  desirable  to  make  a 
measurement  setup  that  will  perform  the  entire  measurement  automatically.  For  this  reason,  a 
measurement  environment  was  set  up  as  depicted  in  Figure  1.  A  central  part  is  the  SR620  Universal 
Time-Interval  Counter  [7].  This  device  plays  two  roles:  on  the  one  hand,  it  is  responsible  for  the 
measurement  of  the  period  of  the  oscillator  under  test,  multiple  periods  respectively.  On  the  other  hand,  it 
is  used  as  a  remote  control  for  the  FPGA  board.  This  board  is  used  as  a  frequency  divider  and  provides 
the  counter  with  a  range  of  frequencies  based  on  the  temperature-stabilized  oscillator  under  test.  In  order 
to  cancel  out  the  influence  of  clock  uncertainties  originating  from  device  internal  oscillators,  the  SR620 
device  was  driven  by  a  stable  rubidium  frequency  standard.  The  overall  measurement  is  controlled  by  a 
PC,  which  also  stores  the  measurement  results. 


Figure  1.  Measurement  setup. 


MEASUREMENT  ACQUISITION 

To  demonstrate  the  described  modeling  approach,  the  Allan  variance  of  a  2  MHz  oscillator  was 
determined.  This  was  done  by  taking  samples  using  different  averaging  times  of  the  oscillator  with  the 
help  of  the  measuring  setup  described  in  the  previous  section.  The  range  of  x  was  chosen  to  be  set  to 
logarithmic  equally  spaced  values  1  ps  to  500  s.  This  upper  value  of  x  is  determined  by  measurement 
equipment.  Now  the  Allan  variance  can  be  calculated  by  [8]: 
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1 

2(iV-l) 


N- 1 


n= 1 


where  yn  is  the  normalized  frequency  departure  defined  as: 

(*„+l  ~Xn) 

y n 

T 


(9) 


(10) 


where  xn  (the  so-called  time  errors)  are  the  measured  samples  and  N  is  the  number  of  samples.  After 
calculating  the  results,  the  Allan  variance  versus  x  is  displayed  on  a  double  logarithmic  plot. 

The  next  step  is  to  calculate  the  values  of  the  h,  parameters.  Although  the  structure  of  the  model  (2) 
suggests  the  straightforward  approach  of  choosing  five  measured  values  and  the  solution  of  the  resulting 
linear  equation  system,  this  turns  out  not  to  be  the  optimal  approach.  The  latter  is  caused  by  the  fact  that, 
in  general,  imperfections  such  as  noise  will  not  lead  to  an  optimal  fit.  That  is  why  an  approximation  of 
the  Allan  variance  and  the  parameters  of  (2)  are  needed.  This  function  needs  to  be  in  form  of  equation 
(4).  The  approximation  can  be  done  in  several  ways,  although  not  all  of  them  have  shown  to  be 
successful  in  the  general  case.  Using  the  methods  of  the  least  quadratic  error  does  not  prove  to  be  very 
successful,  as  the  logarithmic  and,  therefore,  not  equally  spaced  measurements  and  the  relatively  big 
range  of  measurement  values  do  not  lead  to  an  optimal  convergence.  However,  instead  of  using  the  sum 
of  the  quadratic  error,  or  ||  .||2.,  the  infinity  norm  ||  .  H*,  of  the  error  vector  was  used.  However,  the  success 
this  method  heavily  depends  on  the  starting  values  of  this  unconstrained,  nonlinear  approximation;  a  first 
rough  guess  of  these  values  can  be  obtained  using  the  least  quadratic  error  method. 

Figure  2  shows  the  Allan  variance  as  a  function  of  x  displayed  on  a  double  logarithmic  plot  for  two  cases: 
the  red-marked  squares  have  been  determined  by  measurement,  and  the  blue  curve  represents  the 
approximation  using  these  points. 


Fitted  function  and  Measurement  values  of  o2 


Figure  2.  The  Allan  variance  plot  of  the  real  oscillator  and  of  the  approximation. 
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Since  the  approximated  function  has  the  form  of: 

a2y  =  1CT18  •  (0.18974971293808r  +  0.28527328756472  +  0.2271 1027708165^  + 

0. 10828422348094r'2  +  0.00000002073475  r'3)  (11) 

The  h;  parameters  can  be  easily  calculated  with  equation  (7).  The  calculated  values  of  the  h  parameters 
are  illustrated  in  Table  1. 


Table  1.  The  values  of  the  hj  parameters. 


i 

Value 

-2 

1 .2846933726509597  •  10  20 

-1 

2.0580462073483913-10  19 

0 

4.542215104500719-10  19 

1 

5.2216853435436694-10  20 

2 

5.2216853435436694-10  20 

CONCLUSION  AND  FURTHER  WORK 

For  discrete  event  simulations  of  high-precision  clock  synchronization,  a  proper  simulation  model  of 
oscillators  is  more  than  important.  In  fact  a  model,  also  covering  statistical  behavior  such  as  jitter,  and 
not  only  an  aging  frequency  offset,  is  mandatory.  This  turns  out  to  be  not  too  easy,  as  it  has  to  be  possible 
to  recall  the  simulator  model  several  times,  but  no  output  value  must  be  inconsistent  with  a  previous 
answer  of  the  simulation.  Such  a  model  can  be  built  via  filtering  a  white-noise  source  with  a  filter  that 
provides  the  same  Allan  deviation  as  a  real-world  oscillator.  This  paper  explains  how  to  extract  those 
parameters  of  actual  measurements  using  the  infinity  norm  ||  .  f  method.  This  leads  to  the  parameters  of 
the  spectral  power  density  that  can  be  used  as  a  basis  for  the  filter  design.  This  approach  is  demonstrated 
with  some  real-world  measurements  of  an  oscillator. 

As  a  next  step,  this  model  of  the  power  density  spectrum  will  be  estimated  by  a  pink-noise  filter,  and 
more  and  deeper  investigations  concerning  measurement-model  fitting  will  be  done. 
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